Fertilization Analysis (Final)

Load libraries

knitr::opts_chunk$set(echo = TRUE)

list.of.packages <- c("gsheet", "tidyverse", "janitor", "plotly", "glmmTMB", "metafor", "broom.mixed", "car", "ggpubr", "scales") #add new libraries here 

new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)

# Load all libraries 
lapply(list.of.packages, FUN = function(X) {
  do.call("require", list(X)) 
})
sessionInfo()

Prepare data

NOTE: data is read directly from the GoogleSheet using a share link that was set to “anyone with a link can view”

# Read in data from GoogleSheet 
data.fert <- as_tibble(gsheet2tbl('https://docs.google.com/spreadsheets/d/111SuH548Et6HDckjbQjtYk-B8A5VfNkUYZgcUNRPxxY/edit?usp=sharing'))

## Warning: Missing column names filled in: 'X22' [22]

#replace NR (not reported) with NA and convert columns to factor / numeric where needed 
data.fert <- data.fert %>% 
  na_if("NR") %>%     
  clean_names()   %>% # fill in spaces with underscores for column names 
  mutate_at(c('phylum', 'study', 'taxonomic_group', 'common_name', 'latin_name', 'error_statistic'), as.factor) %>%
  mutate_at(c('p_h_experimental', 'p_h_control', 'p_co2_experimental', 'p_co2_control', 'ave_fert_percent_p_h', 'error_percent_p_h', 'number_trials_p_h', 'insemination_time', 'sperm_per_m_l', 'sperm_egg_ratio', 'number_females', 'number_males'), as.numeric)

## Warning: NAs introduced by coercion

## Warning: NAs introduced by coercion

## Warning: NAs introduced by coercion

## Warning: NAs introduced by coercion

## Warning: NAs introduced by coercion

#data.fert %>% select(study, ave_fert_percent_p_h, error_percent_p_h, error_statistic, number_trials_p_h) %>% View()

Explore data with figures

ggplotly(
data.fert %>%
ggplot(mapping=aes(x=p_h_experimental, y=ave_fert_percent_p_h, group=taxonomic_group, col=taxonomic_group, text=`common_name`)) + 
  geom_point(size=1.5, width=0.02) +
  facet_wrap(~phylum, scale="free") +
  geom_smooth(method="lm", se=TRUE, aes(fill=taxonomic_group)) +
  ggtitle("Fertilization Rate ~ pH exposure by phylum, linear") +
  theme_minimal() 
)

## Warning: Ignoring unknown parameters: width

## `geom_smooth()` using formula 'y ~ x'

## Warning: Removed 25 rows containing non-finite values (stat_smooth).

ggplotly(
data.fert %>%
ggplot(mapping=aes(x=p_h_experimental, y=ave_fert_percent_p_h, group=study, col=study, text=`common_name`)) + 
  geom_point(size=1.5, width=0.02) +
  facet_wrap(~phylum, scale="free") +
  geom_smooth(method="lm", se=TRUE, aes(fill=taxonomic_group)) +
  ggtitle("Fertilization Rate ~ pH exposure by phylum, linear") +
  theme_minimal() 
)

## Warning: Ignoring unknown parameters: width

## `geom_smooth()` using formula 'y ~ x'

## Warning: Removed 25 rows containing non-finite values (stat_smooth).

## Warning in qt((1 - level)/2, df): NaNs produced

## Warning in qt((1 - level)/2, df): NaNs produced

## Warning in qt((1 - level)/2, df): NaNs produced

## Warning in qt((1 - level)/2, df): NaNs produced

## Warning in qt((1 - level)/2, df): NaNs produced

## Warning in qt((1 - level)/2, df): NaNs produced

## Warning in qt((1 - level)/2, df): NaNs produced

## Warning in qt((1 - level)/2, df): NaNs produced

ggplotly(
data.fert %>%
ggplot(mapping=aes(x=p_h_experimental, y=ave_fert_percent_p_h, group=taxonomic_group, col=taxonomic_group, text=`common_name`)) + 
  geom_point(size=1, width=0.02) +
  facet_wrap(~phylum, scale="free") +
  geom_smooth(method="lm", se=TRUE, formula=y ~ poly(x, 2, raw=TRUE), aes(fill=taxonomic_group)) +
  ggtitle("Fertilization Rate ~ pH exposure by phylum, polynomial") +
  theme_minimal()
)

## Warning: Ignoring unknown parameters: width

## Warning: Removed 25 rows containing non-finite values (stat_smooth).

## Warning in predict.lm(model, newdata = new_data_frame(list(x = xseq)), se.fit =
## se, : prediction from a rank-deficient fit may be misleading

Convert all error values to separate SE & SD columns

95% Confidence Interval

Upper 95% CI = Mean + 1.96*SE; I recorded the difference between the upper 95%CI and the mean (Upper 95%CI - Mean). To convert I will use:

SE = (Upper 95%CI - Mean) / 1.96
SD = ((Upper 95%CI - Mean) / 1.96) * sqrt(n)

Standard Deviation / Standard Error conversions

SE= SD/sqrt(n)
SD = SE*sqrt(n)
where n=sample size

data.fert <- data.fert %>% 
  mutate(SE =  case_when(error_statistic == "SD" ~ error_percent_p_h/sqrt(number_trials_p_h), 
         error_statistic == "95% CI" ~ error_percent_p_h/1.96,
         is.na(error_statistic) ~ error_percent_p_h, 
         error_statistic == "SE" ~ error_percent_p_h)) 

data.fert <- data.fert %>% 
  mutate(SD =  case_when(error_statistic == "SE" ~ error_percent_p_h*sqrt(number_trials_p_h), 
         error_statistic == "95% CI" ~ (error_percent_p_h/1.96)*sqrt(number_trials_p_h),
         is.na(error_statistic) ~ error_percent_p_h, 
         error_statistic == "SD" ~ error_percent_p_h)) 

Calculate delta pH (pH experimental - pH control)

data.fert <- data.fert %>% 
  mutate(pH_delta = p_h_control-p_h_experimental)

Convert % fertilization data to proportion fertilized, and replace any values >1 (aka 100% fertilized) with 1

data.fert <- data.fert %>% 
  mutate(ave_fert_proport = case_when(ave_fert_percent_p_h <= 100 ~ ave_fert_percent_p_h/100,
                                      ave_fert_percent_p_h > 100 ~ 1))

Inspect data

data.fert$ave_fert_proport %>% hist()

ggplot(data.fert, aes(group=phylum, col=phylum)) + geom_density(aes(ave_fert_proport))

## Warning: Removed 13 rows containing non-finite values (stat_density).

# How many studies per ~phylum? 
data.fert %>%
  select(phylum, study) %>%
  distinct(phylum, study) %>%
  group_by(phylum) %>% count()

## # A tibble: 4 x 2
## # Groups:   phylum [4]
##   phylum            n
##   <fct>         <int>
## 1 Cnidaria          6
## 2 Crustacea         5
## 3 Echinodermata    12
## 4 Mollusca         18

# How many studies per taxonomic group? 
data.fert %>%
  select(phylum, study, taxonomic_group) %>%
  distinct(phylum, study, taxonomic_group) %>%
  group_by(taxonomic_group) %>% count()

## # A tibble: 10 x 2
## # Groups:   taxonomic_group [10]
##    taxonomic_group     n
##    <fct>           <int>
##  1 abalone             2
##  2 clam                5
##  3 copepod             3
##  4 coral               6
##  5 mussel              4
##  6 non-copepod         2
##  7 non-urchin          2
##  8 oyster              5
##  9 scallop             3
## 10 sea urchin         11

Calculate weights for models

weights <- metafor::escalc(measure='MN',
                mi=data.fert$ave_fert_percent_p_h,
                sdi = data.fert$SD,
                ni=data.fert$number_trials_p_h, options(na.action="na.pass"))  
data.fert$w <-weights$vi

How many studies per phylum after filtering out those w/o weight calculation?

# How many studies per ~phylum? 
data.fert %>% drop_na(w) %>%
  select(phylum, study) %>%
  distinct(phylum, study) %>%
  group_by(phylum) %>% count()

## # A tibble: 4 x 2
## # Groups:   phylum [4]
##   phylum            n
##   <fct>         <int>
## 1 Cnidaria          5
## 2 Crustacea         5
## 3 Echinodermata    11
## 4 Mollusca         15

Generate binomial models

• Include study as random effect for all models
  
• Include varying combinations of phylum & experimental pH
  
• Decided NOT to include taxa, since so few studies per group

Test all candidate covariates alone - are any significant? Answer: only experimental pH

Anova(glmmTMB(ave_fert_proport ~ p_h_experimental + (1|study), data=drop_na(data.fert, w), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #yes

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: ave_fert_proport
##                   Chisq Df Pr(>Chisq)   
## p_h_experimental 9.1418  1   0.002498 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Anova(glmmTMB(ave_fert_proport ~ p_h_control + (1|study), data=drop_na(data.fert, w), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #no

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: ave_fert_proport
##              Chisq Df Pr(>Chisq)
## p_h_control 0.0811  1     0.7758

Anova(glmmTMB(ave_fert_proport ~ phylum + (1|study), data=drop_na(data.fert, w), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #no

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: ave_fert_proport
##         Chisq Df Pr(>Chisq)
## phylum 1.9538  3      0.582

Anova(glmmTMB(ave_fert_proport ~ insemination_time + (1|study), data=drop_na(data.fert, w), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #no

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: ave_fert_proport
##                    Chisq Df Pr(>Chisq)
## insemination_time 0.0557  1     0.8134

Anova(glmmTMB(ave_fert_proport ~ sperm_per_m_l + (1|study), data=drop_na(data.fert, w), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #no

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: ave_fert_proport
##               Chisq Df Pr(>Chisq)
## sperm_per_m_l 0.024  1      0.877

Anova(glmmTMB(ave_fert_proport ~ sperm_egg_ratio + (1|study), data=drop_na(data.fert, w), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #no

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): Model convergence problem; non-positive-
## definite Hessian matrix. See vignette('troubleshooting')

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: ave_fert_proport
##                 Chisq Df Pr(>Chisq)
## sperm_egg_ratio 2e-04  1       0.99

Anova(glmmTMB(ave_fert_proport ~ number_females + (1|study), data=drop_na(data.fert, w), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #no

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: ave_fert_proport
##                Chisq Df Pr(>Chisq)
## number_females 0.054  1     0.8162

Anova(glmmTMB(ave_fert_proport ~ number_males + (1|study), data=drop_na(data.fert, w), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #no

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: ave_fert_proport
##               Chisq Df Pr(>Chisq)
## number_males 0.2089  1     0.6476

Test all candidate covariates with experimental pH - are any significant? Answer: only experimental pH

Anova(glmmTMB(ave_fert_proport ~ p_h_experimental*p_h_control + (1|study), data=drop_na(data.fert, w), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #no

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: ave_fert_proport
##                               Chisq Df Pr(>Chisq)   
## p_h_experimental             9.1034  1   0.002551 **
## p_h_control                  0.0432  1   0.835310   
## p_h_experimental:p_h_control 0.0804  1   0.776795   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Anova(glmmTMB(ave_fert_proport ~ p_h_experimental*insemination_time + (1|study), data=drop_na(data.fert, w), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #no

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: ave_fert_proport
##                                      Chisq Df Pr(>Chisq)   
## p_h_experimental                   10.6640  1   0.001092 **
## insemination_time                   0.7586  1   0.383757   
## p_h_experimental:insemination_time  0.8087  1   0.368510   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Anova(glmmTMB(ave_fert_proport ~ p_h_experimental*phylum + (1|study), data=drop_na(data.fert, w), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #no

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: ave_fert_proport
##                          Chisq Df Pr(>Chisq)   
## p_h_experimental        7.9631  1   0.004774 **
## phylum                  4.2472  3   0.235975   
## p_h_experimental:phylum 2.9060  3   0.406354   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Anova(glmmTMB(ave_fert_proport ~ p_h_experimental*sperm_per_m_l + (1|study), data=drop_na(data.fert, w), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #no

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): Model convergence problem; non-positive-
## definite Hessian matrix. See vignette('troubleshooting')

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: ave_fert_proport
##                                 Chisq Df Pr(>Chisq)
## p_h_experimental               1.1881  1     0.2757
## sperm_per_m_l                  0.0272  1     0.8691
## p_h_experimental:sperm_per_m_l 0.0009  1     0.9762

Anova(glmmTMB(ave_fert_proport ~ p_h_experimental*sperm_egg_ratio + (1|study), data=drop_na(data.fert, w), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #no

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): Model convergence problem; non-positive-
## definite Hessian matrix. See vignette('troubleshooting')

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: ave_fert_proport
##                                   Chisq Df Pr(>Chisq)  
## p_h_experimental                 3.6744  1    0.05526 .
## sperm_egg_ratio                  0.0047  1    0.94507  
## p_h_experimental:sperm_egg_ratio 0.0000  1    0.99590  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Anova(glmmTMB(ave_fert_proport ~ p_h_experimental*number_females + (1|study), data=drop_na(data.fert, w), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #no

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: ave_fert_proport
##                                  Chisq Df Pr(>Chisq)   
## p_h_experimental                8.7577  1   0.003083 **
## number_females                  0.0718  1   0.788736   
## p_h_experimental:number_females 0.0786  1   0.779138   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Anova(glmmTMB(ave_fert_proport ~ p_h_experimental*number_males + (1|study), data=drop_na(data.fert, w), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #no

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: ave_fert_proport
##                                Chisq Df Pr(>Chisq)   
## p_h_experimental              8.8671  1   0.002903 **
## number_males                  0.1893  1   0.663490   
## p_h_experimental:number_males 0.0179  1   0.893581   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Save best fit model (called “best”) to object, and examine (experimental pH as sole predictor)

car::Anova(best <- glmmTMB(ave_fert_proport ~ p_h_experimental + (1|study), data=drop_na(data.fert, w), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #phylum, pH, & phylum:pH sign. factors 

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: ave_fert_proport
##                   Chisq Df Pr(>Chisq)   
## p_h_experimental 9.1418  1   0.002498 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(best)

##  Family: binomial  ( logit )
## Formula:          ave_fert_proport ~ p_h_experimental + (1 | study)
## Data: drop_na(data.fert, w)
## Weights: 1/(w + 1)^2
## 
##      AIC      BIC   logLik deviance df.resid 
##     73.7     84.8    -33.9     67.7      295 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  study  (Intercept) 1.929    1.389   
## Number of obs: 298, groups:  study, 33
## 
## Conditional model:
##                  Estimate Std. Error z value Pr(>|z|)   
## (Intercept)       -22.382      7.673  -2.917  0.00354 **
## p_h_experimental    3.098      1.025   3.023  0.00250 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Write a function to predict %fertilization at various pH levels (and specifying the model)

predict.fert <- function(pH, model) {
    linear.predictor <- as.vector((summary(best)$coefficients$cond[,"Estimate"]["(Intercept)"] + 
                           summary(best)$coefficients$cond[,"Estimate"]["p_h_experimental"]*pH))
    predicted <- exp(linear.predictor) / (1+exp(linear.predictor))
    return(paste("Fertilization rate predicted for pH ", pH, ": ", 
                 scales::percent(x=predicted, accuracy = .01), sep=""))
} 

# Estimate % fertilization @ pH 8.0 using hand-typed equation 
exp(-22.382104 + 3.098070*8) / (1 + exp(-22.382104 + 3.098070*8)) 

## [1] 0.9170144

# now use function 
predict.fert(8.0, best)

## [1] "Fertilization rate predicted for pH 8: 91.70%"

predict.fert(7.5, best)

## [1] "Fertilization rate predicted for pH 7.5: 70.13%"

predict.fert(7.0, best)

## [1] "Fertilization rate predicted for pH 7: 33.28%"

predict.fert(6.0, best)

## [1] "Fertilization rate predicted for pH 6: 2.20%"

Generate estimates & confidence intervals (log likelihood)

confint(best)

##                                      2.5 %    97.5 %   Estimate
## cond.(Intercept)               -37.4217617 -7.342447 -22.382104
## cond.p_h_experimental            1.0897904  5.106349   3.098070
## study.cond.Std.Dev.(Intercept)   0.5899624  3.270399   1.389033

Inspect residuals ~ fitted values

aa5 <- augment(best, data=drop_na(data.fert, w))

## Warning in mr - fitted(object): longer object length is not a multiple of
## shorter object length

#ggplotly(
ggplot(aa5, aes(x=.fitted,y=.resid)) + 
    geom_point() + 
  geom_smooth()

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

## Warning: Removed 12 rows containing non-finite values (stat_smooth).

## Warning: Removed 12 rows containing missing values (geom_point).

#)

Generate model predictions and plot against real data

ph.min.max <- drop_na(data.fert, w) %>% 
  select(phylum, p_h_experimental) %>% 
  group_by(phylum) %>% 
  summarize(min=min(p_h_experimental, na.rm=TRUE), max=max(p_h_experimental, na.rm=TRUE))
  
phylum.list <- list()
for (i in 1:nrow(ph.min.max)) {
  phylum.list[[i]] <- data.frame(ph=c(seq(from=as.numeric(ph.min.max[i,"min"]), 
                            to=as.numeric(ph.min.max[i,"max"]), 
                            by=0.01)),
                   phylum=rep(c(ph.min.max[i,"phylum"])))
}
new.data <- bind_rows(phylum.list) %>% purrr::set_names(c("p_h_experimental", "phylum"))
new.data$study <- NA 
new.data$w <- NA 

predict.test.df <- predict(best, newdata = new.data, se.fit = TRUE, type="response")
predict.test.df.df <- predict.test.df %>%
  as.data.frame() %>%
  cbind(new.data)

#scales::show_col(c("#e41a1c","#4daf4a","#ff7f00","#984ea3",'#377eb8'))

Figure caption:

Fertilization success (%) by experimental pH across marine taxa examined in this review. Meta-analysis was performed using a binomial regression model, and indicates that fertilization success decreases with pH across Crustacea (5 studies), Echinodermata (12 studies), and Mollusca (18 studies). Fertilization success was not significantly affected by pH in Cnidaria (4 studies). Each point reflects the average % fertilization reported by one study at an experimental pH.

Generate single figure with all phyla

# Examine pH-experimental significance as predictor using X-squared test  
Anova(best)

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: ave_fert_proport
##                   Chisq Df Pr(>Chisq)   
## p_h_experimental 9.1418  1   0.002498 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#ggplotly(
ggplot() + 
  geom_jitter(data=drop_na(data.fert, w), aes(x=p_h_experimental, y=ave_fert_proport, group=phylum, col=phylum, label=study), size=1.2, width=0.03) + #, col="gray40"
  #facet_wrap(~phylum, scales="free") + 
  theme_minimal() +
  ggtitle("Fertilization success ~ pH with binomial-regression model predictions (all phyla)") + 
  xlab("Experimental pH") + ylab("Proportion fertilization success") +
  scale_color_manual(name=NULL, values=c("#ca0020","#f4a582","#92c5de",'#0571b0')) + 
  theme_minimal() + 
  coord_cartesian(ylim = c(0, 1), xlim =c(6,8.5)) +
  geom_line(data = predict.test.df.df, aes(x=p_h_experimental, y=fit), col="gray50") + #, col=phylum
  geom_ribbon(data = predict.test.df.df, aes(x=p_h_experimental, ymin=fit-se.fit, ymax=fit+se.fit), linetype=2, alpha=0.1, col="gray50") + theme(legend.position = "top") + guides(colour = guide_legend(override.aes = list(size=4))) +
  annotate(geom="text", x=6.4, y=0.9, size=5, colour="gray50", label=paste("χ2 p-value =", round(Anova(best)[,"Pr(>Chisq)"], digits =4), sep=" ")) #) #add fill=phylum in geom_ribbon aes if color desired 

## Warning: Ignoring unknown aesthetics: label

## Warning: Removed 12 rows containing missing values (geom_point).

Examine color scheme (show_col is from scales library)

show_col(c("#ca0020","#f4a582","#92c5de",'#0571b0'))

#e41a1c = Cnidaria (red)  
#ff7f00 = Crustacea (orange)  
#4daf4a = Echinodermata (green)  
#377eb8 = Mollusca (blue)  

#ca0020 = coral 
#f4a582 = crustacean 
#92c5de = echinoderm
#0571b0 = mollusc

Generate figure: pH data by phylum, each with binomial regression model fit + CI

ggplotly(
ggplot() + 
  geom_jitter(data=drop_na(data.fert, w), aes(x=p_h_experimental, y=ave_fert_proport, group=phylum, col=phylum, label=study), size=1.2, width=0.03, col="gray40") +
  facet_wrap(~phylum, scales="free") + theme_minimal() +
  ggtitle("Fertilization success ~ pH with binomial-regression model predictions") + 
  xlab("Experimental pH") + ylab("Proportion fertilization success") +
  #scale_color_manual(name=NULL, values=c("#e41a1c","#ff7f00","#4daf4a",'#377eb8')) +
  theme_minimal() + 
  coord_cartesian(ylim = c(0, 1), xlim =c(6,8.5)) +
  geom_line(data = predict.test.df.df, aes(x=p_h_experimental, y=fit), col="gray50") + #, col=phylum
  geom_ribbon(data = predict.test.df.df, aes(x=p_h_experimental, ymin=fit-se.fit, ymax=fit+se.fit), linetype=2, alpha=0.1, col="gray50") + theme(legend.position = "none")) #add fill=phylum in geom_ribbon aes if color desired 

## Warning: Ignoring unknown aesthetics: label

Run GLMs on each phylum and generate plots

1. Are slopes significantly different from zero?
2. If so, what is the equation?

Cnidaria

Anova(glmmTMB(ave_fert_proport ~ p_h_experimental + (1|study), data=data.fert %>% drop_na(w) %>% filter(phylum=="Cnidaria"), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #no

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: ave_fert_proport
##                   Chisq Df Pr(>Chisq)
## p_h_experimental 1.5152  1     0.2183

Anova(glmmTMB(ave_fert_proport ~ p_h_control + (1|study), data=data.fert %>% drop_na(w) %>% filter(phylum=="Cnidaria"), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #no

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: ave_fert_proport
##             Chisq Df Pr(>Chisq)
## p_h_control 1.555  1     0.2124

Anova(glmmTMB(ave_fert_proport ~ insemination_time + (1|study), data=data.fert %>% drop_na(w) %>% filter(phylum=="Cnidaria"), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #no

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: ave_fert_proport
##                    Chisq Df Pr(>Chisq)
## insemination_time 0.1092  1     0.7411

Anova(glmmTMB(ave_fert_proport ~ sperm_per_m_l + (1|study), data=data.fert %>% drop_na(w) %>% filter(phylum=="Cnidaria"), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #no

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): Model convergence problem; non-positive-
## definite Hessian matrix. See vignette('troubleshooting')

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: ave_fert_proport
##               Chisq Df Pr(>Chisq)
## sperm_per_m_l        1

Anova(glmmTMB(ave_fert_proport ~ sperm_egg_ratio + (1|study), data=data.fert %>% drop_na(w) %>% filter(phylum=="Cnidaria"), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #no

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): Model convergence problem; non-positive-
## definite Hessian matrix. See vignette('troubleshooting')

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: ave_fert_proport
##                 Chisq Df Pr(>Chisq)
## sperm_egg_ratio        1

Anova(glmmTMB(ave_fert_proport ~ number_females + (1|study), data=data.fert %>% drop_na(w) %>% filter(phylum=="Cnidaria"), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #no

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: ave_fert_proport
##                Chisq Df Pr(>Chisq)
## number_females 4e-04  1     0.9844

Anova(glmmTMB(ave_fert_proport ~ number_males + (1|study), data=data.fert %>% drop_na(w) %>% filter(phylum=="Cnidaria"), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #no

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: ave_fert_proport
##               Chisq Df Pr(>Chisq)
## number_males 0.0058  1     0.9392

# None significant, but still develop model for plot 
model.cnidaria <- glmmTMB(ave_fert_proport ~ p_h_experimental + (1|study), data=data.fert %>% drop_na(w) %>% filter(phylum=="Cnidaria"), binomial(link = "logit"), na.action=na.exclude,  weights = 1/(w+1)^2)

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

Anova(model.cnidaria)

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: ave_fert_proport
##                   Chisq Df Pr(>Chisq)
## p_h_experimental 1.5152  1     0.2183

summary(model.cnidaria) 

##  Family: binomial  ( logit )
## Formula:          ave_fert_proport ~ p_h_experimental + (1 | study)
## Data: data.fert %>% drop_na(w) %>% filter(phylum == "Cnidaria")
## Weights: 1/(w + 1)^2
## 
##      AIC      BIC   logLik deviance df.resid 
##     28.1     34.0    -11.1     22.1       50 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  study  (Intercept) 0.4241   0.6513  
## Number of obs: 53, groups:  study, 5
## 
## Conditional model:
##                  Estimate Std. Error z value Pr(>|z|)
## (Intercept)       -41.697     34.526  -1.208    0.227
## p_h_experimental    5.406      4.392   1.231    0.218

paste("Experimental pH χ2 p-value =", round(Anova(model.cnidaria)[,"Pr(>Chisq)"], digits = 3), sep=" ")

## [1] "Experimental pH χ2 p-value = 0.218"

predict.cnidaria <- predict(model.cnidaria, newdata = subset(new.data, phylum=="Cnidaria"), se.fit = TRUE, type="response")
predict.cnidaria.df <- predict.cnidaria %>%
  as.data.frame() %>%
  cbind(subset(new.data, phylum=="Cnidaria"))

plot.cnidaria <- ggplot() + 
  geom_jitter(data=data.fert %>% drop_na(w) %>% 
                filter(phylum=="Cnidaria"), aes(x=p_h_experimental, y=ave_fert_proport), size=1.2, width=0.03, col="#ca0020") +
  ggtitle("Cnidaria") + 
  xlab(NULL) + ylab("Proportion fertilization success") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5, size=13, colour="gray30")) +
  coord_cartesian(ylim = c(0, 1), xlim =c(6,8.5)) +
  geom_line(data = predict.cnidaria.df, aes(x=p_h_experimental, y=fit), col="#ca0020") + 
  geom_ribbon(data = predict.cnidaria.df, aes(x=p_h_experimental, ymin=fit-se.fit, ymax=fit+se.fit), linetype=2, alpha=0.1, col="#ca0020") + 
  annotate(geom="text", x=6.5, y=0.85, size=3.5, colour="gray20", label=paste("χ2 p-value =\n", round(Anova(model.cnidaria)[,"Pr(>Chisq)"], digits = 3), sep=" ")) #Significance of experimental pH \nas predictor: 
print(plot.cnidaria)

Crustacea

Anova(glmmTMB(ave_fert_proport ~ p_h_experimental + (1|study), data=data.fert %>% drop_na(w) %>% filter(phylum=="Crustacea"), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #no

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in (function (start, objective, gradient = NULL, hessian = NULL, : NA/
## NaN function evaluation

## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: ave_fert_proport
##                   Chisq Df Pr(>Chisq)
## p_h_experimental 0.0635  1      0.801

Anova(glmmTMB(ave_fert_proport ~ p_h_control + (1|study), data=data.fert %>% drop_na(w) %>% filter(phylum=="Crustacea"), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #no

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in (function (start, objective, gradient = NULL, hessian = NULL, : NA/
## NaN function evaluation

## Warning in (function (start, objective, gradient = NULL, hessian = NULL, : NA/
## NaN function evaluation

## Warning in (function (start, objective, gradient = NULL, hessian = NULL, : NA/
## NaN function evaluation

## Warning in (function (start, objective, gradient = NULL, hessian = NULL, : NA/
## NaN function evaluation

## Warning in (function (start, objective, gradient = NULL, hessian = NULL, : NA/
## NaN function evaluation

## Warning in (function (start, objective, gradient = NULL, hessian = NULL, : NA/
## NaN function evaluation

## Warning in (function (start, objective, gradient = NULL, hessian = NULL, : NA/
## NaN function evaluation

## Warning in (function (start, objective, gradient = NULL, hessian = NULL, : NA/
## NaN function evaluation

## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')

## Warning in fitTMB(TMBStruc): Model convergence problem; function evaluation
## limit reached without convergence (9). See vignette('troubleshooting')

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: ave_fert_proport
##              Chisq Df Pr(>Chisq)
## p_h_control 0.0028  1     0.9579

Anova(glmmTMB(ave_fert_proport ~ insemination_time + (1|study), data=data.fert %>% drop_na(w) %>% filter(phylum=="Crustacea"), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #no

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: ave_fert_proport
##                   Chisq Df Pr(>Chisq)
## insemination_time        1

model.crustacea <- glmmTMB(ave_fert_proport ~ p_h_experimental + (1|study), data=data.fert %>% drop_na(w) %>% filter(phylum=="Crustacea"), binomial(link = "logit"), na.action=na.exclude,  weights = 1/(w+1)^2)

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in (function (start, objective, gradient = NULL, hessian = NULL, : NA/
## NaN function evaluation

## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')

Anova(model.crustacea)

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: ave_fert_proport
##                   Chisq Df Pr(>Chisq)
## p_h_experimental 0.0635  1      0.801

summary(model.crustacea) 

##  Family: binomial  ( logit )
## Formula:          ave_fert_proport ~ p_h_experimental + (1 | study)
## Data: data.fert %>% drop_na(w) %>% filter(phylum == "Crustacea")
## Weights: 1/(w + 1)^2
## 
##      AIC      BIC   logLik deviance df.resid 
##       NA       NA       NA       NA       21 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance   Std.Dev.
##  study  (Intercept) 7.022e-109 8.38e-55
## Number of obs: 24, groups:  study, 5
## 
## Conditional model:
##                  Estimate Std. Error z value Pr(>|z|)
## (Intercept)        -20.51      75.22  -0.273    0.785
## p_h_experimental     2.58      10.24   0.252    0.801

predict.crustacea <- predict(model.crustacea, newdata = subset(new.data, phylum=="Crustacea"), se.fit = TRUE, type="response")
predict.crustacea.df <- predict.crustacea %>%
  as.data.frame() %>%
  cbind(subset(new.data, phylum=="Crustacea"))

plot.crustacea <- ggplot() + 
  geom_jitter(data=data.fert %>% drop_na(w) %>% 
                filter(phylum=="Crustacea"), aes(x=p_h_experimental, y=ave_fert_proport), size=1.2, width=0.03, col="#f4a582") +
  ggtitle("Crustacea") + 
  xlab(NULL) + ylab(NULL) +
  theme_minimal() + 
  theme(plot.title = element_text(hjust = 0.5, size=13, colour="gray30")) +
  coord_cartesian(ylim = c(0, 1), xlim =c(6,8.5)) +
  geom_line(data = predict.crustacea.df, aes(x=p_h_experimental, y=fit), col="#f4a582") +
  geom_ribbon(data = predict.crustacea.df, aes(x=p_h_experimental, ymin=fit-se.fit, ymax=fit+se.fit), linetype=2, alpha=0.1, col="#f4a582")  + 
  annotate(geom="text", x=6.5, y=0.85, size=3.5, colour="gray20", label=paste("χ2 p-value =\n", round(Anova(model.crustacea)[,"Pr(>Chisq)"], digits = 3), sep=" "))
print(plot.crustacea)

Echinodermata

Anova(glmmTMB(ave_fert_proport ~ p_h_experimental + (1|study), data=data.fert %>% drop_na(w) %>% filter(phylum=="Echinodermata"), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #yes

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: ave_fert_proport
##                   Chisq Df Pr(>Chisq)  
## p_h_experimental 6.3225  1    0.01192 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Anova(glmmTMB(ave_fert_proport ~ p_h_control + (1|study), data=data.fert %>% drop_na(w) %>% filter(phylum=="Echinodermata"), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #no

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: ave_fert_proport
##              Chisq Df Pr(>Chisq)
## p_h_control 0.2848  1     0.5936

Anova(glmmTMB(ave_fert_proport ~ insemination_time + (1|study), data=data.fert %>% drop_na(w) %>% filter(phylum=="Echinodermata"), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #no

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: ave_fert_proport
##                    Chisq Df Pr(>Chisq)  
## insemination_time 3.7605  1    0.05248 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Anova(glmmTMB(ave_fert_proport ~ sperm_per_m_l + (1|study), data=data.fert %>% drop_na(w) %>% filter(phylum=="Echinodermata"), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #no

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: ave_fert_proport
##               Chisq Df Pr(>Chisq)
## sperm_per_m_l     0  1     0.9963

Anova(glmmTMB(ave_fert_proport ~ sperm_egg_ratio + (1|study), data=data.fert %>% drop_na(w) %>% filter(phylum=="Echinodermata"), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #no

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: ave_fert_proport
##                  Chisq Df Pr(>Chisq)
## sperm_egg_ratio 0.0671  1     0.7957

Anova(glmmTMB(ave_fert_proport ~ number_females + (1|study), data=data.fert %>% drop_na(w) %>% filter(phylum=="Echinodermata"), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #no

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: ave_fert_proport
##                 Chisq Df Pr(>Chisq)
## number_females 0.2046  1      0.651

Anova(glmmTMB(ave_fert_proport ~ number_males + (1|study), data=data.fert %>% drop_na(w) %>% filter(phylum=="Echinodermata"), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #no

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: ave_fert_proport
##               Chisq Df Pr(>Chisq)
## number_males 2.4125  1     0.1204

model.echinodermata <- glmmTMB(ave_fert_proport ~ p_h_experimental + (1|study), data=data.fert %>% drop_na(w) %>% filter(phylum=="Echinodermata"), binomial(link = "logit"), na.action=na.exclude,  weights = 1/(w+1)^2)

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

Anova(model.echinodermata)

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: ave_fert_proport
##                   Chisq Df Pr(>Chisq)  
## p_h_experimental 6.3225  1    0.01192 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(model.echinodermata) 

##  Family: binomial  ( logit )
## Formula:          ave_fert_proport ~ p_h_experimental + (1 | study)
## Data: data.fert %>% drop_na(w) %>% filter(phylum == "Echinodermata")
## Weights: 1/(w + 1)^2
## 
##      AIC      BIC   logLik deviance df.resid 
##     29.0     37.9    -11.5     23.0      142 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  study  (Intercept) 2.661    1.631   
## Number of obs: 145, groups:  study, 11
## 
## Conditional model:
##                  Estimate Std. Error z value Pr(>|z|)  
## (Intercept)       -41.938     16.894  -2.482   0.0131 *
## p_h_experimental    5.619      2.235   2.514   0.0119 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

predict.echinodermata <- predict(model.echinodermata, newdata = subset(new.data, phylum=="Echinodermata"), se.fit = TRUE, type="response")
predict.echinodermata.df <- predict.echinodermata %>%
  as.data.frame() %>%
  cbind(subset(new.data, phylum=="Echinodermata"))

plot.echino <- ggplot() + 
  geom_jitter(data=data.fert %>% drop_na(w) %>% 
                filter(phylum=="Echinodermata"), aes(x=p_h_experimental, y=ave_fert_proport), size=1.2, width=0.03, col="#92c5de") +
  ggtitle("Echinodermata") + 
  xlab("Experimental pH") + ylab("Proportion fertilization success") +
  theme_minimal() + 
  theme(plot.title = element_text(hjust = 0.5, size=13, colour="gray30")) +
  coord_cartesian(ylim = c(0, 1), xlim =c(6,8.5)) +
  geom_line(data = predict.echinodermata.df, aes(x=p_h_experimental, y=fit), col="#92c5de") +
  geom_ribbon(data = predict.echinodermata.df, aes(x=p_h_experimental, ymin=fit-se.fit, ymax=fit+se.fit), linetype=2, alpha=0.1, col="#92c5de") + 
  annotate(geom="text", x=6.4, y=0.85, size=3.5, colour="gray20", label=paste("χ2 p-value =\n", round(Anova(model.echinodermata)[,"Pr(>Chisq)"], digits = 3), sep=" "))
print(plot.echino)

Mollusca

Anova(glmmTMB(ave_fert_proport ~ p_h_experimental + (1|study), data=data.fert %>% drop_na(w) %>% filter(phylum=="Mollusca"), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #meh 

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: ave_fert_proport
##                   Chisq Df Pr(>Chisq)  
## p_h_experimental 2.8595  1    0.09084 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Anova(glmmTMB(ave_fert_proport ~ p_h_control + (1|study), data=data.fert %>% drop_na(w) %>% filter(phylum=="Mollusca"), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #no

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in (function (start, objective, gradient = NULL, hessian = NULL, : NA/
## NaN function evaluation

## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: ave_fert_proport
##             Chisq Df Pr(>Chisq)
## p_h_control        1

Anova(glmmTMB(ave_fert_proport ~ insemination_time + (1|study), data=data.fert %>% drop_na(w) %>% filter(phylum=="Mollusca"), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #no

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: ave_fert_proport
##                   Chisq Df Pr(>Chisq)
## insemination_time 0.008  1     0.9288

Anova(glmmTMB(ave_fert_proport ~ sperm_per_m_l + (1|study), data=data.fert %>% drop_na(w) %>% filter(phylum=="Mollusca"), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #no

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): Model convergence problem; non-positive-
## definite Hessian matrix. See vignette('troubleshooting')

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: ave_fert_proport
##               Chisq Df Pr(>Chisq)
## sperm_per_m_l        1

Anova(glmmTMB(ave_fert_proport ~ sperm_egg_ratio + (1|study), data=data.fert %>% drop_na(w) %>% filter(phylum=="Mollusca"), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #no

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in eval(family$initialize): Model convergence problem; non-positive-
## definite Hessian matrix. See vignette('troubleshooting')

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: ave_fert_proport
##                 Chisq Df Pr(>Chisq)
## sperm_egg_ratio        1

Anova(glmmTMB(ave_fert_proport ~ number_females + (1|study), data=data.fert %>% drop_na(w) %>% filter(phylum=="Mollusca"), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #no

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: ave_fert_proport
##                 Chisq Df Pr(>Chisq)
## number_females 0.5883  1     0.4431

Anova(glmmTMB(ave_fert_proport ~ number_males + (1|study), data=data.fert %>% drop_na(w) %>% filter(phylum=="Mollusca"), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #no

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: ave_fert_proport
##               Chisq Df Pr(>Chisq)
## number_males 0.3482  1     0.5552

model.mollusca <- glmmTMB(ave_fert_proport ~ p_h_experimental + (1|study), data=data.fert %>% drop_na(w) %>% filter(phylum=="Mollusca"), binomial(link = "logit"), na.action=na.exclude,  weights = 1/(w+1)^2)

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

Anova(model.mollusca)

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: ave_fert_proport
##                   Chisq Df Pr(>Chisq)  
## p_h_experimental 2.8595  1    0.09084 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(model.mollusca) 

##  Family: binomial  ( logit )
## Formula:          ave_fert_proport ~ p_h_experimental + (1 | study)
## Data: data.fert %>% drop_na(w) %>% filter(phylum == "Mollusca")
## Weights: 1/(w + 1)^2
## 
##      AIC      BIC   logLik deviance df.resid 
##     19.8     26.8     -6.9     13.8       73 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance  Std.Dev.
##  study  (Intercept) 2.209e-06 0.001486
## Number of obs: 76, groups:  study, 14
## 
## Conditional model:
##                  Estimate Std. Error z value Pr(>|z|)  
## (Intercept)       -9.5492     6.4052  -1.491   0.1360  
## p_h_experimental   1.5196     0.8986   1.691   0.0908 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

predict.mollusca <- predict(model.mollusca, newdata = subset(new.data, phylum=="Mollusca"), se.fit = TRUE, type="response")
predict.mollusca.df <- predict.mollusca %>%
  as.data.frame() %>%
  cbind(subset(new.data, phylum=="Mollusca"))

plot.mollusca <- ggplot() + 
  geom_jitter(data=data.fert %>% drop_na(w) %>% 
                filter(phylum=="Mollusca"), aes(x=p_h_experimental, y=ave_fert_proport), size=1.2, width=0.03, col="#0571b0") +
  ggtitle("Mollusca") + 
  xlab("Experimental pH") + ylab(NULL) +
  theme_minimal() + 
  theme(plot.title = element_text(hjust = 0.5, size=13, colour="gray30")) +
  coord_cartesian(ylim = c(0, 1), xlim =c(6,8.5)) +
  geom_line(data = predict.mollusca.df, aes(x=p_h_experimental, y=fit), col="#0571b0") + 
  geom_ribbon(data = predict.mollusca.df, aes(x=p_h_experimental, ymin=fit-se.fit, ymax=fit+se.fit), linetype=2, alpha=0.1, col="#0571b0")  + 
  annotate(geom="text", x=6.4, y=0.9, size=3.5, colour="gray20", label=paste("χ2 p-value =\n", round(Anova(model.mollusca)[,"Pr(>Chisq)"], digits = 3), sep=" "))
print(plot.mollusca)

## Warning: Removed 12 rows containing missing values (geom_point).

all.plots <- ggarrange(plot.cnidaria + ylab(NULL), plot.crustacea, plot.echino + ylab(NULL), plot.mollusca + rremove("y.text"), ncol=2, nrow=2)

## Warning: Removed 12 rows containing missing values (geom_point).

annotate_figure(all.plots, left = text_grob("Proportion Fertilization Success", color = "gray20", rot = 90))

#add the following to remove x and y axis labels: ` + rremove("x.text") + rremove("y.text")`
# add the following to add plot labels: `, labels=c("A", "B", "C", "D")`

Additional analysis steps, for posterity: inspect candidate covariates - are these important experimental design factors?

convert candidate covariates to numeric

data.fert <- data.fert %>%
  mutate_at(vars(number_trials_p_h, insemination_time, sperm_per_m_l, sperm_egg_ratio, number_females, number_males), as.numeric) 

Test Control pH

car::Anova(covariates1 <- glmmTMB(ave_fert_proport ~ p_h_control*p_h_experimental, data=drop_na(data.fert, w), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #no interaction 

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: ave_fert_proport
##                               Chisq Df Pr(>Chisq)   
## p_h_control                  0.0001  1   0.992776   
## p_h_experimental             9.4009  1   0.002169 **
## p_h_control:p_h_experimental 0.7993  1   0.371315   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ggplotly(
  ggplot(data.fert, aes(x=p_h_control, y=ave_fert_proport)) + 
    geom_jitter(aes(color = p_h_experimental, label=common_name), size=2.2, pch=19) + ggtitle("% Fertilization ~ Control pH\nsign. interaction") + 
    scale_color_gradientn(colours = rainbow(5))) 

## Warning: Ignoring unknown aesthetics: label

#    scale_color_gradient(low = "red", high = "blue"))

Test Insemination time

car::Anova(covariates2 <- glmmTMB(ave_fert_proport ~  insemination_time*p_h_experimental, data=drop_na(data.fert, w), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #no interaction  

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: ave_fert_proport
##                                      Chisq Df Pr(>Chisq)    
## insemination_time                   0.7835  1  0.3760640    
## p_h_experimental                   10.8791  1  0.0009725 ***
## insemination_time:p_h_experimental  0.7018  1  0.4021760    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ggplotly(
  ggplot(data.fert, aes(x=insemination_time, y=ave_fert_proport)) + 
    geom_point(aes(color = p_h_experimental, label=common_name), size=2.2, pch=19) + ggtitle("% Fertilization ~ Insemination Time\nsign. interaction_") + 
    scale_color_gradientn(colours = rainbow(5))) 

## Warning: Ignoring unknown aesthetics: label

#    scale_color_gradient(low = "red", high = "blue"))

Test sperm concentration

car::Anova(covariates3 <- glmmTMB(ave_fert_proport ~ sperm_per_m_l*p_h_experimental, data=drop_na(data.fert, w), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #not sign 

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: ave_fert_proport
##                                 Chisq Df Pr(>Chisq)
## sperm_per_m_l                  0.0121  1     0.9124
## p_h_experimental               1.3803  1     0.2400
## sperm_per_m_l:p_h_experimental 0.0004  1     0.9849

ggplotly(ggplot(data.fert, aes(x=sperm_per_m_l, y=ave_fert_proport)) + 
           geom_point(aes(color = p_h_experimental, label=common_name), size=2.2, pch=19) + 
           ggtitle("% Fertilization ~ Sperm Concentration\nnot sign.") + 
          scale_color_gradientn(colours = rainbow(5))) 

## Warning: Ignoring unknown aesthetics: label

#           scale_color_gradient(low = "red", high = "blue"))

Test sperm:egg ratio

car::Anova(covariates4 <- glmmTMB(ave_fert_proport ~ sperm_egg_ratio*p_h_experimental, data=drop_na(data.fert, w), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #not sign 

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: ave_fert_proport
##                                   Chisq Df Pr(>Chisq)   
## sperm_egg_ratio                  0.0073  1   0.931826   
## p_h_experimental                 6.8239  1   0.008994 **
## sperm_egg_ratio:p_h_experimental 0.0000  1   0.995566   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ggplotly(ggplot(data.fert, aes(x=sperm_egg_ratio, y=ave_fert_proport)) + 
           geom_point(aes(color = p_h_experimental, label=common_name), size=2.2, pch=19) + 
           ggtitle("% Fertilization ~ Sperm:Egg Ratio\nnot sign.") + 
    scale_color_gradientn(colours = rainbow(5))) 

## Warning: Ignoring unknown aesthetics: label

#           scale_color_gradient(low = "red", high = "blue"))

Test no. females used for experimental eggs

car::Anova(covariates5 <- glmmTMB(ave_fert_proport ~  number_females*p_h_experimental, data=drop_na(data.fert, w), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #not sign.

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: ave_fert_proport
##                                  Chisq Df Pr(>Chisq)   
## number_females                  0.0001  1   0.991430   
## p_h_experimental                9.2586  1   0.002344 **
## number_females:p_h_experimental 0.0550  1   0.814625   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ggplotly(ggplot(data.fert, aes(x=number_females, y=ave_fert_proport)) + 
           geom_point(aes(color = p_h_experimental, label=common_name), size=2.2, pch=19) + 
           ggtitle("% Fertilization ~ No. Females\nno interaction") + 
    scale_color_gradientn(colours = rainbow(5))) 

## Warning: Ignoring unknown aesthetics: label

#           scale_color_gradient(low = "red", high = "blue"))

Test no. males used for experimental eggs

car::Anova(covariates6 <- glmmTMB(ave_fert_proport ~  number_males*p_h_experimental, data=drop_na(data.fert, w), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #not sign. 

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: ave_fert_proport
##                                Chisq Df Pr(>Chisq)   
## number_males                  0.0535  1   0.817025   
## p_h_experimental              9.2455  1   0.002361 **
## number_males:p_h_experimental 0.0423  1   0.836968   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ggplotly(ggplot(data.fert, aes(x=number_males, y=ave_fert_proport)) + 
           geom_point(aes(color = p_h_experimental, label=common_name), size=2.2, pch=19) + 
           ggtitle("% Fertilization ~ No. Males\nsign. interaction") + 
     scale_color_gradientn(colours = rainbow(5))) 

## Warning: Ignoring unknown aesthetics: label

#          scale_color_gradient(low = "red", high = "blue"))